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Highly excited states for isospin | baryons are calculated for the first time using lattice QCD 
with two flavors of dynamical quarks. Anisotropic lattices are used with two pion masses, m n = 
416(36) MeV and 578(29) MeV. The lowest four energies are reported in each of the six irreducible 
representations of the octahedral group at each pion mass. The lattices used have dimensions 
24 3 x 64, spatial lattice spacing a B « 0.11 fm and temporal lattice spacing a t = ^a s . Clear evidence 

is found for a | state in the pattern of negative-parity excited states. This agrees with the pattern 
of physical states and spin | has been realized for the first time on the lattice. 

PACS numbers: 11.15.Ha,12.38.Gc,12.38.Lg 



I. INTRODUCTION 

A major goal for lattice quantum chromodynamics (QCD) is the determination of the spectrum of the excited 
states of QCD. This goal is a complement to experimental work that studies the hadrons and their excitations and 
decays. In recent years, large amounts of data have been collected at Jefferson Laboratory regarding the spectrum of 
excitations of nucleons. The Excited Baryon Analysis Center aims to analyze the data using the best hadronic models 
available. [H 0] Lattice QCD calculations are needed as a means to link this program to the Lagragian of QCD. 

When QCD was formulated as the basic theory that would explain hadrons and their excited states, it could not 
be solved for the mass spectrum from first principles because of its nonperturbative nature. Much effort over the past 
thirty years has been devoted to developing the methods and tools to solve QCD on a lattice. Accurate resolution 
of the excited states of hadrons using lattice QCD has proven difficult. In Euclidean space, excited state correlation 
functions decay faster than the ground state. At large times, the signals for excited states are swamped by the signals 
for lower energy states. Improved resolution in the temporal direction is essential for progress. An anisotropic lattice 
where the temporal lattice spacing is finer than spatial spacings can provide better resolution while avoiding the 
increase in computational cost associated with a similar reduction of all spacings. The improved resolution must be 
combined with two other ingredients. A large number of operators is required that overlap well with excited states. 
The use of variational methods is essential to separate the excited states. 

Large sets of baryon operators were developed and projected to the irreducible representations of the octahedral 
group in Refs. [HQ Link smearing and quark smearing were found to both be needed in order to optimize the quality 
of the signals that are obtained with the operators in Ref. [j| . Variational methods were used to determine the spectra 
of I = h and I = | excited baryons using the quenched approximation in Ref. [|| . 

In this work, we take another step toward the goal of determining the spectrum of nucleon excited states by studying 
the spectrum of isospin i excited baryons in two-flavor QCD, using u and d quarks that have the same mass. Results 
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are obtained on 24 3 x64 lattices with two values of the pion mass: 416(36) MeV and 578(29) MeV. 

In the physical spectrum for isospin i the lowest three states are the nucleon, N, the Roper resonance, N' (Pn) 
and the opposite-parity N* (Six)). Quenched lattice QCD calculations 0, B H, E3, El, E3] generally have found a 
spectrum inverted with respect to experiment, with the N' heavier than the negative-parity N*. An exception is 
the calculation of the Kentucky group [l3j| that obtained the correct mass ordering with a pion mass below 400 MeV 
(after subtracting the effects of the quenched "ghosts"). This has helped to motivate full QCD simulations where 
the spectrum can be determined without unphysical contributions from "ghost" states. Moreover, many additional 
excited states have been observed experimentally that should be reproduced by lattice QCD and full-QCD simulations 
are needed as a complement to the experimental searches for new excited states. 

Anisotropic techniques have been adopted in lattice calculations for relativistic heavy quark actions for the spectrum 
of charmonium Q3, [l5| , for calculations of the spectrum of glueballs [l6| and to extract excited baryon states [1, 
EL H [13, EH, EH, HH ■ Previous results using anisotropic lattices include two- flavor anisotropic dynamical simulations 
performed by CP-PACS [H[ and the TrinLat collaboration [II]. 

In this work, we report the nucleon spectrum using the interpolating basis from Refs. on two-flavor, anisotropic, 
Wilson fermion and Wilson gauge configurations. The action parameters and bare gauge and fermion anisotropies 
are tuned such that the gauge anisotropy (as determined from Wilson loop ratios) and the fermion anisotropy (as 
determined from the meson dispersion relation) are both consistent with the desired renormalizcd anisotropy a s /a t = 3. 
Our configurations were generated using the Chroma (23[ HMC code with multi-timescale integration. 

The organization of this paper is as follows: In Sec.m we discuss the details of the actions used and their parameters. 
Then in Sec. lIIII we discuss the Hybrid Monte Carlo (HMC) used in this work and show how it is applied to anisotropic 
lattices with mass preconditioning. Section IIVI presents results for the conventional determination of hadron masses 
and the anisotropy from two-point correlation functions and Sec. |V] presents our procedure and results for setting 
the lattice scale in physical units. Section [VII discusses the construction of large numbers of baryon operators in the 
relevant irreducible representations of the octahedral group and demonstrates the noise suppression that is obtained 
by smearing both the quark and gauge fields. Section IVIII presents results for the 1 = \ baryon spectrum for pion 

masses of 416 MeV and 578 MeV using Nf = 2 lattices. Clear evidence for a spin-| state is presented. Some 
conclusions are presented in Sec. IVIIII 



II. LATTICE ACTIONS 



In this section, we describe the gauge and fermion actions used in this calculation. For the gauge sector, we use a 
Wilson anisotropic gauge action 



sk[U] 



ft 



n PBsr ( x ) + J2$n-pA* 



where flw = ReTr(l — V) and V is the plaquette 

P„ v (x) = U^(x)U„(x + rfUUx + u)Ut{x). 



(1) 



(2) 



The coupling g 2 appears in j3 = 2N c /g 2 . The parameter £o is the bare gauge anisotropy. In the fermion sector, wc 
adopt the anisotropic Wilson fermion action [2J] 



Sf[[/,V,V>] - a 3 s a t ^2^{x)M w ^(x), 

X 

M w = mo + v t W t + u s W a , 



(3) 



where 



Wu = 



V„/(x) 
A,f(x) 



2a L 



UpWfix + ri-Ukx-rifix-n) 
U,(x)f(x + /i) + UUx - fx)f(x - /i) - 2f(x) 



(4) 
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In terms of dimensionless variables = a^ 2 ip, mo = moa t , V M = a^V^, A M = a^A M , and the dimensionless "Wilson 
operator" = — g-y^A^, we find that the fermion matrix My/ becomes 

M w = —\ a t mo + VtWt + £ Yl W.\ ■ (5) 

Because it is possible to redefine the fields as in Refs. [25|,[26|, one coefficient (either v t or v s ) is redundant; here we set 
Vt = 1 and v s — v for tuning. For convenience of parameterization, we use the bare gauge and fermion anisotropics, 
7 S) /, defined as 

7 9 = £o, 7/ = "^- (6) 

The parameters 7 ff , 7^ and the quark mass mo require tuning in order to realize the desired renormalization 
constraints. The bare gauge and fermion anisotropy parameters 7 g and 7/ are tuned to obtain the desired renormalized 
gauge and fermion anisotropics (£ s and both equal to a s /a t = 3.0. The renormalized gauge anisotropy (£ s ) can 
be determined from the static-quark potential using Klassen's "Wilson- loop ratio" [27[ : 

iy ss (a; + l,y) 



(7) 



where W s t are the Wilson loops involving the temporal direction, and W ss are those involving only the spatial 
directions. We determine the renormalized gauge anisotropy £ g by minimizing pll ] 

r/n _r {Rss{x,y) - R st {x^ g y)) 2 

where Ai? s and Ai? t are the statistical errors of R ss and R s t. A fixed background gauge field in the spatial "z" 
direction is used following the Schrodinger-functional scheme p8| which allows for a determination of the critical mass 
using the PCAC Ward identity. For more details, see Sec. IV B of Ref. [29(. We determine the renormalized fermion 
anisotropy £/ through the conventional relativistic meson dispersion relation as will be discussed in Sec. IIV Al 

We find that when £ ff = £ = 2.38, u = 1 (or = £ = 2.38), £ g and £/ (see Sec. IIV A[) are consistent with 3, given 
our other choices of parameters. The critical mass at these bare parameters is m c = —0.41473. The mo parameter 
within our range of interest has negligible effect on the anisotropics. (Similar results are observed in the three-flavor 
anisotropic clover action study in Ref. [2{|.) We set mo to —0.4086 and —0.4125 for our pion-mass study. 



III. ALGORITHM 

The anisotropic Wilson configurations were generated with the Hybrid Monte Carlo (HMC) algorithm [3(J. To in- 
crease the efficiency of the method we employed several techniques such as Hasenbusch st yle mass preconditioning [3ll ] , 
the use of multiple timescale integration schemes [Hf , chronological inversion methods 1331 during the molecular dy- 
namics, and evolving the temporal links with different time-steps than the spatial ones [22|]. We discuss some of the 
pertinent details below: 



A. Hybrid Monte Carlo 

The basic technique for gauge generation is a Markov Chain Monte Carlo method where one moves from an 
initial gauge configuration to a successive one by generating a new trial configuration and then performing an accep- 
tance/rejection test upon the new one. If the trial configuration is accepted, it becomes the next configuration in the 
chain, otherwise the original configuration becomes the next state in the chain. 

In order to use a global Metropolis accept/reject step with a reasonable acceptance rate, the space of states is 
extended to include momenta tt^{x) canonical to the gauge links U^(x) so that one may define a Hamiltonian 



(9) 
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where S is the action. It is then possible to propose new configurations from previous ones by performing Hamiltonian 
Molecular Dynamics (MD) to get from the initial to the proposed state. Using a reversible and area preserving MD 
evolution maintains detailed balance, which is sufficient for the algorithm to converge. In order to ensure ergodicity 
in the entire phase space, the momenta need to change periodically. This can be accomplished by refreshing the 
momenta from a Gaussian heat bath prior to the MD update step. 

In order to deal with the fermion determinant, it is standard to use the method of pseudofermions. One integrates 
out the Grassman- valued fermion fields in the action and rewrites the resulting determinant as an integral over bosonic 
fields, 

Z = J [df 1 ][dr 1 ]e- ,1 ' D ' 1 = det (V) = J [d# t ][d0e-* tl '~ 1 * J (10) 

where rj and fj are the Grassman valued fields, T> is some Hcrmitian, positive-definite kernel and 4^ and 4> are the 
bosonic pseudofermion fields. Our phase space is thus enlarged to include also the pseudofermion fields. Like the 
momenta, these fields need to be refreshed before each MD step to carry out the pseudofermion integral. 
In the case of a two-flavor simulation, T> is typically of the form, 

V = Q^Q. (11) 

For the rest of this work Q is an even-odd preconditioned fermion matrix for an individual flavor of fermion. In this 
case T> is manifestly Hermitian and positive definite, and the integral in Eq. [10] is guaranteed to exist. Furthermore, 
the pseudofermion fields can easily be refreshed by producing a vector \ filled with Gaussian noise with a variance of 
i and then forming <f> = Q^X- 



B. Multiple Time Scale Anisotropic Molecular Dynamics Update 

While any reversible and area-preserving MD update scheme can be used in the MD step, the acceptance rate is 
controlled by the truncation error in the scheme. This manifests itself as a change in the Hamiltonian, SH, over an 
MD trajectory, since we use the Metropolis acceptance probability 

P acc = mm(l,e- SH ) . (12) 

We may easily construct a reversible scheme by combining symplectic update steps U p {8t) and U q (Sr) which update 
momenta and coordinates by a time step of length 5r respectively 

Uri&rp): Mx),U M (x)) -> (ir^x) + F^x)6t^U^(x)) , (13) 
U q {8^): (7r M (x),l7 M (x)) -> (tt m (s), e^^U^x)) , (14) 

where F ll (x) is the MD force coming from the variation of the action with respect to the gauge fields. We emphasize 
that one may update all the links pointing in direction /i with a separate step size St^. While this may not be useful 
in isotropic simulations, in an anisotropic calculation with one fine direction, it may be advantageous to use a shorter 
timcstep to update the links in that direction to ameliorate the typically larger forces that result from the shorter 
lattice spacing [22l | . The anisotropy in step size requires a small amount of manual fine tuning, but should be similar 
to the anisotropy in the lattice spacings. 

Our base integration scheme in this work is due to Omelyan [H, H3, [35| ; we use the combined update operator 

^(St) = U p {XST)U q (^ST)U p (l - 2X6T)U q (^5T)Up(X6T), (15) 

which results in a scheme that is clearly reversible and is accurate to 0(5t 3 ). The size of the leading error term can 
be further minimized by tuning the parameter A. In our work we used the value of A from Ref. [34| without any 
further tuning, which promises an efficiency increase of approximately 50% over the simple leapfrog algorithm. 

In Refs. [33, Hi| it was shown that a reversible, multi-level integration scheme can be constructed which allows 
various pieces of the Hamiltonian to be integrated at different timescales. Let us consider a Hamiltonian of the form 



H(n, U) = -7Tl(x)n,(x) + S^U) + S 2 (U), 



(16) 
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where S\(U) and 52(f) are pieces of the action with corresponding MD forces F\ and F 2 respectively. One can then 
split the integration into 2 timescales. One can integrate with respect to action Si(U) using U 1 (5ti), where in the 
component U p (Srx) we use only the force Fx. The whole system can then be integrated with the update 

U 2 (St 2 ) =U' p {\8t 2 )U 1 {] 1 5t 2 )U' p {\ - 2\6t 2 )U 1 (^5t 2 )u; } (\5t 2 ), (17) 

where in U p we update the momenta using only F 2 . Thus we end up with two characteristic integration timescales 
Stx and St 2 . The scheme generalizes recursively to a larger number of scales. A criterion for tuning the algorithm 
is to arrange for terms in the action to be mapped to different timescales so that on two timescales i and j we have 
||-Fi||<fa} ~ ||-Fj||<5T,, as suggested in Rcf. [3l| . We now proceed to outline how we split our action. 
We can write our gauge action schematically as 

S = S S (U) + S t (U), (18) 

where the term 5 S contains only loops with spatial gauge links, and the 5 t term contains loops with spatial and 
temporal links. While the term 5 S produces forces only in the spatial directions, the St term produces forces in both 
the spatial and the temporal directions. In particular the spatial forces from St are larger in magnitude than the 
spatial forces from 5 S by roughly the order of the anisotropy, and in turn, the temporal forces from St arc larger 
than the spatial forces from St- Our anisotropic integration step size balances the spatial and temporal forces of the 
St term against each other. However, in order to balance the spatial forces from St and 5 S against each other, we 
integrate them on separate time scales. 



C. Mass Preconditioning 

Following the work of [3l| our fermion determinant for the two flavor simulation can be written as 

. + , det (QtQ) . . 

det (QtQ) = p, *\ det (Q{Q h ) , (19) 

det [QlQt 



where Q is the fermion matrix with our desired fermion mass m and Qh is the fermion matrix for which we choose a 
heavier fermion mass to/j. After introducing pseudofermions the fermion action can be written as 

Sf =5} +Sj, (20) 

where 

5} = (j)\Qh (Q^Q)' 1 Qlfa, (21) 

S 2 f = 4 (QlQh) ~V (22) 

This trick introduces two main advantages: first, because Qh is heavier than Q, inversion in S'j will take fewer 
iterations than solving with Q directly, and forces resulting from Qh will likewise be smaller than those that would 
result from Q allowing slightly longer time steps; second, as long as to is not very different from rrih, we have that 
to first order Qh (Q^Q) Q\ ~ 1 + A and that fluctuations with gauge fields will be to first order given by It 
should be clear, that as rrih — > to we have A — * 0, and that the resulting force F — ► 0, in other words, that the 
magnitude of the force from Sj can be made small in a controlled manner. The result is that while the inversions in 
Sj are performed with Q and can be quite costly; by choosing rrih appropriately the force from Sj can be reduced so 
that Sj can be put on a long time scale and evaluated relatively infrequently during an MD trajectory. Some amount 
of effort is required to tune rrih so that the number of force evaluations from Sj can be suitably reduced, while at 
the same time keeping rrih heavy enough, so that the force evaluations and inversions from 5 2 do not become overly 
expensive. 



D. Chronological Inversion Methods 



In order to further reduce our inversion costs, we employ chronological guesses in our MD. Before every new solve 
we produce a chronological guess by employing the Minimal Residual Extrapolation method (MRE) of [HI]. This 
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method works by using the last n solution vectors, which are orthonormalizcd with respect to each other to create an 
?i dimensional basis. Let us denote these basis vectors u,-. The new initial guess v g is then constructed as 

v g = 2J ami i = l...n, (23) 

i 

where aj are coefficients to be determined given the constraint that the resulting v g minimize the functional minimized 
by the Conjugate Gradients process in the subspace spanned by tj^: 

* M = vlQlQv g - xS - v\ X , (24) 

where \ IS the right-hand side of the linear system for which the initial guess is being generated. Minimizing the 
functional ^ with respect to <2j leads to the following set of linear equations for af. 

N 
i=l 

We emphasize that the use of chronological solution methods introduces reversibility violations into the MD evolu- 
tion, and so the equations must be solved essentially exactly to avoid the reversibility violations from becoming large, 
and affecting the detailed balance condition and thereby biasing the Monte Carlo Markov process. To this end in our 
simulations we required a relative stopping residuum of 

rMD . IH-WOMI < 10 -». (26) 
llxll 



E. Summary 



In summary, our HMC algorithm uses a Hamiltonian composed of the kinetic piece, the two gauge action pieces S 3 
and St and the two fermion action pieces and Sj. Our Molecular Dynamics evolution uses a 2nd Order Omelyan 
integrator split over three timescales: 

• Time scale 1 is the slowest, with time step 8t\, and is used to evolve the Hasenbusch ratio term with action S^; 

• Time scale 2 is faster, with time step 8t%, and it is used to evolve the mass preconditioned fermion term with 
action S'j and the spatial gauge term with action S s ; 

• Time scale 3 is the fastest, with time step 8t$, and it is used to evolve the temporal gauge term with action St. 

Our overall MD trajectory length is set to be r = 1.0. In addition at all levels of the integrator, the spatial time step 
on that level is a factor of £,md = 2.4 larger than the temporal step. Both fermionic terms use the MRE chronological 
guess technique with up to n — 8 previous solutions. These preconditioning masses rrih and the concrete step sizes 
are summarized in Table HI 

The acceptance rates were typically between 60% and 70%. The simulations at mass m = —0.4125 made use of the 
QCDOC supercomputer [l4j, as well as BlueGene Teragrid Resource at San Diego Supercomputer Center, while the 
entire m = —0.4086 dataset was generated on Jaguar, a Cray XT3 resource at the National Center for Computational 
Science (NCCS) at Oak Ridge National Laboratory through the INCITE'07 program The HMC algorithm with the 
various imp rovements discussed in this section is implemented and is freely available as part of the Chroma software 
system |23j . 



IV. CONVENTIONAL SPECTROSCOPY 



A. Meson spectrum 



We use meson interpolating fields of the form i/jTifi and i/jFj^i/j that overlap with the physical states listed in 
Table ITT] Correlation functions are calculated and we fit them with the analytical function, 

C(t) = A[e- mt + e -™( T -*)) , (27) 
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p 


mo 


mh 


5ti 


Sr 2 
8r-\ 


<St 3 
5ra 


5.5 


-0.4086 


-0.3700 


i 

4 


1 

2 


1 

3 


5.5 


-0.4125 


-0.3740 


1 

4 


1 
4 


1 

3 



TABLE I: The mass preconditioning masses mj, and the time steps used in the Omelyan integration scheme in our simulations, 
for each target sea quark mass mo- Except for time scale 1, the time step for each time scale is given relative to the previous 
one. Trajectories are of length r = 1, with a step size anisotropy of (,md = 2.4. The target solver residuum was r — 10~ 8 for 
both MD and Energy calculations and each fermionic term employed the MRE chronological guess method with up to the last 
8 previous vectors. 
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I = 1 


0"+ 


75 


7T 


1~ 


7m 


P 


0++ 


1 


fflo 


1++ 


7m 75 


ai 


1+- 


7m 7^ 


bi 



TABLE II: Meson interpolating operators. The indicated charge-conjugation (C) quantum numbers apply only to particles 
with zero net flavor. 

where m is the mass parameter and T is the time extent of the lattice. Results for the mass parameters obtained from 
the fits of the correlations functions are summarized in Tables IIIII and IIVI Comparisons of the fits for the 7r meson 
with the effective mass are shown in Figure [T] for the case of quark mass parameter moa t = -0.4125. Horizontal lines 
show the corresponding pion mass parameter of the fits and the error band. 

The fermion anisotropy is determined through the conventional relativistic meson dispersion relation: 

E 2 (p) = m 2 + (28) 

where the energy E(p) and the mass m are in units of at, and p = ^p 1 where L s is the spatial lattice size in units 
of a s . From the two-point correlation functions we calculate the energy E at the spatial momenta p = ^p 1 for 
n = (0, 0, 0), (1, 0, 0), (1, 1, 0) and (2, 0, 0) (averaged). The fitted jackknife energies are used in a linear fit of E 2 (p) as 
a function of p 2 as in Eq. (|2"8|) in order to extract £/. Figure shows the dispersion relations for ir and p mesons. The 
fitted values of £/ are 2.979(28) (with first three momenta fit) for the ir meson and 3.045(35) (with first four momenta 
fit) for the rho meson. The central (green) line shows the fit and the blue bands show the errors. The desired fermion 
anisotropy matches the gauge anisotropy, which is 3 in our case. 
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Sink 
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Pi 


0.0754(5) 




S 


Si 


0.0748(5) 


0.1430(7) 


s 


P 2 


0.0747(7) 


0.1438(10) 


s 


s 2 


0.0747(5) 


0.1427(8) 


s 


P1&P2 


0.0753(5) 


0.1431(7) 


s 


S1&S2 


0.0747(5) 


0.1427(8) 


Result 


0.0750(7) 


0.1431(8) 



TABLE III: Meson masses (in temporal lattice units) for Nf — 2 with light quark mass moat = —.4125 based on 862 
configurations. Columns 1 and 2 label the sources and sinks as smeared (S) or point (P) and the correlation functions based 
on the operators of Table ITT1 (Si or Pi), or based on including an extra factor 7* in the operators of Table HT1 (S2 or P2). 
Simultaneous fits are performed for two types of operators in the results of rows 5 and 6. Fits are performed in time windows 
of (26-32) with x 2 /DOF = 0.45(15) for the tt meson. 
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Source 


Sink 




Trip 


S 


Pi 


0.1088(7) 


0.1668(12) 


S 


Si 




0.1652(13) 


S 


P 2 


0.1088(12) 


0.1680(13) 


S 


P1&S1 




0.1670(11) 


S 


P1&P2 


0.1088(7) 




Results 


0.1088(8) 


0.1668(16) 



TABLE IV: Meson masses (in temporal lattice units) for Nf = 2 with quark mass moa t = —.4086 based on 363 configurations. 
Notation is the same as in TableHni Fits are performed in time windows of (26-32) with x 2 /DOF = 0.45(15) for the ■K meson. 



0.1 



0.095 - 




5 10 15 20 25 30 35 

t 

FIG. 1: Pion effective mass results for different smearings (S denotes smeared, P denotes point) and operators (opl denotes 
operators of Table Ull and op2 denotes operators including an extra factor of 7*). Results are based on 862 gauge configurations, 
quark mass parameter mo = -0.4125 and a 24 3 x 64 lattice. 



V. SCALE SETTING 



In order to set the scale, the heavy-quark (static) potential V(r) is calculated on a 16 3 x64 lattice. This is expected 
to have the form 

V(r) =C +- +crr, (29) 
r 




u 1 1 1 1 1 1 1 1 

0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 
P 2 

FIG. 2: Dispersion relation for tt (left panel) and p (right panel) 
lattice. 
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P 2 

for quark mass parameter mo = -0.4125 and 24 3 x64 
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where r is the separation between the static quarks. The scale implied by the heavy quark potential is often specified 
using the Sommer parameter r which is defined by the condition 



,dV(r) 



dr 



1.65. 



On the lattice, we calculate Wilson loops which determine the static quark potential via 

W(r,t) = Ae~ v(r)t . 



(30) 



(31) 



In order to improve the signal and to extract the potential V(r) from smaller time separations, we smear the gauge 
links in the spatial directions using stout smearing with parameters n p = 16, n p p = 2.5. We fit the Wilson loops as a 
function of t at each available r to determine V(r). We further fit V(r) to determine C, a and a according to Eq. (|29|) 
using a standard jackknife procedure. Putting these parameters back into Eq. (|30j) . we solve for rn/a«. 

Finally, we relate this to the physical scale by using the value ro= 0.462(H)(4) fm from Refs. [33, [38[ and set the 
scale a,. The results are summarized in Table W\ 



ro(fm) 


miat 


ro/cis 


a s (fm) 


aT/^MeV) 




m^MeV) 


Co 


0.462(H)(4) 
0.462(H)(4) 


-0.4086 
-0.4125 


4.10(8) 
4.26(12) 


0.113(7) 
0.108(7) 


5310(265) 
5556(333) 


0.1088(37) 
0.0750(24) 


578(29) 
416(36) 


2.38 
2.38 



TABLE V: The value of the Sommer parameter, ro, is listed in column 1 and the ratio ro/a s for each quark mass on our 16 3 x64 
lattices is listed in column 3. The scale obtained from ro/(ro/a s ) is listed in column 4. Using the renormalized anisotropy £ = 
3, we find the temporal spacing a^ 1 as given in column 5. The pion mass in lattice units is given in column 6 and in MeV units 
in column 7, while the bare anisotropy £o is given in column 8. 



VI. BARYON OPERATORS 



The use of operators whose temporal correlation functions attain their asymptotic form as quickly as possible is 
crucial for reliably extracting excited hadron masses. An important ingredient in constructing such hadron operators 
is the use of smeared fields. Operators constructed from smeared fields have dramatically reduced mixings with the 
high frequency modes of the theory. Both link-smearing and quark-field smearing are necessary. Since excited hadrons 
are expected to be large objects, the use of spatially extended operators is another key ingredient in the operator 
design and implementation. 



A. Smearing 



Spatial links can be smeared using the stout-link procedure described in Ref. (39j . The stout-link smearing scheme 
is analytic, efficient, and produces smeared links that automatically are elements of SU(3) without the need for a 
projection back into SU(3). Note that only spatial staples are used in the link smoothing; no temporal staples are 
used, and the temporal link variables are not smeared. The smeared quark fields can be defined by 

#z)= + \{x), (32) 

where a s and n a are tunable parameters (n a is a positive integer) and the three-dimensional covariant Laplacian 
operators are defined in terms of the smeared link variables Uj(x) as follows: 

AO(x) = (u k (x)0(x + k) - 0(x)\ (33) 

fc=±l,±2,±3^ ' 

where 0{x) is an operator defined at lattice site x with appropriate color structure, and noting that V r -k(%) = Ul(x-k). 

The smeared fields ip and ip arc Grassmann-valucd; in particular, these fields anticommute in the same way that the 
original fields do, and the square of each smeared field vanishes. 
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single- singly- doubly- doubly- triply- triply- 

site displaced displaced-I displaced-L displaced-T displaced-O 

FIG. 3: The spatial arrangements of the extended three-quark baryon operators. Smeared quark-fields are shown by solid 
circles, line segments indicate gauge-covariant displacements, and each hollow circle indicates the location of a Levi-Civita color 
coupling. For simplicity, all displacements have the same length in an operator. Results presented here used displacement 
lengths of 3a s (~ 0.3 fm). 



B. Group theory 

Hadron states are identified by their momentum p, intrinsic spin J, projection A of this spin onto some axis, parity 
P = ±1, and quark flavor content (isospin, strangeness, etc.). Some mesons also include G-parity as an identifying 
quantum number. If one is interested only in the masses of these states, one can restrict attention to the p = sector, 
so operators must be invariant under all spatial translations allowed on a cubic lattice. The little group of all symmetry 
transformations on a cubic lattice which leave p = invariant is the octahedral point group Oh , so operators may be 
classified using the irreducible representations (irreps) of Oh- For mesons, there are ten irreducible representations 
A\ g , A2g, Eg,T\ g ,T2g, A\u, A2U1 E u ,Ti u ,Ti u . The representations with a subscript g(u) are even (odd) under parity. 
The A irreps are one dimensional, the E irreps are two dimensional, and the T irreps are three-dimensional. The 
A\ irreps contain the J = 0, 4, 6, 8, . . . states, the A2 irreps contain the J = 3, 6, 7, 9, . . . states, the E irreps contain 
the J = 2, 4, 5, 6, 7, . . . states, the T\ irreps contain the spin J = 1, 3, 4, 5, . . . mesons, and the T2 irreps contain the 
spin J = 2,3,4,5,... states. For baryons, there are four two-dimensional irreps Gi g ,Gi u ,G2 g , G2 U an d two four- 
dimensional representations H g and H u . The G\ irrep contains the J = i, |, | , . . . states, the H irrep contains 
the J = |, |, |, |, . . . states, and the G2 irrep contains the J = |, |, . . . states. The continuum-limit spins J of 
our states must be deduced by examining degeneracy patterns across the different Oh irreps. 

C. Operator construction and pruning 

Our operators are constructed in a three-stage approach^. First, basic building blocks are chosen. These are taken 
to be smeared covariantly-displaced quark fields 

Q&f")L> -3<J<3, (34) 

where A is a flavor index, a is a color index, a is a Dirac spin index, and the p-link gauge-covariant displacement 
operator in the j-th direction is defined by 

D<f ) (x,x') = U j (x) U 3 (x +])... U^x+ip-l)]^,^, D { p \x,x') = S xxl , (35) 

for j = ±1,±2,±3 and p > 1, and where j ' = defines a zero-displacement operator to indicate no displacement. 
Next, elemental operators Bf(t,x) are devised having the appropriate flavor structure characterized by isospin, 
strangeness, etc., and color structure constrained by gauge invariance. For zero momentum states, translational 
mvanance is imposed: B[(t) = T, x B f(^ x )- Finally, group-theoretical projections arc applied to obtain operators 
which transform irrcducibly under all lattice rotation and reflection symmetries: 

gAXFQ = _^a_ T^(R) U R Bf{t) Ul (36) 



c - u u w 

single- singly- doubly- triply- triply- 

site displaced displaced-L displaced-U displaced-O 

FIG. 4: The spatial arrangements of the quark-antiquark meson operators. In the illustrations, the smeared quarks fields are 
depicted by solid circles, each hollow circle indicates a smeared "barred" antiquark field, and the solid line segments indicate 
covariant displacements. 
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FIG. 5: Effective masses M(i) for unsmeared (black circles) and smeared (red triangles) operators Oss, Osd, Otdt, which are 
representative single-site, singly-displaced, and triply-displaced-T nucleon operators, respectively. Top row: only quark-field 
smearing n a = 32, a s — 4.0 is used. Middle row: only link-variable smearing n p — 16, n p p = 2.5 is applied. Bottom row: 
both quark and link smearing n a = 32, a a = 4.0, n p = 16, n p p — 2.5 are used, dramatically improving the signal for all three 
operators. Results are based on 50 quenched configurations on a 12 3 x 48 anisotropic lattice using the Wilson action with 
a a ~ 0.1 fm, a s /a t ~ 3.0. 



where Off is the double group of Oh, R denotes an element of Off , g r> is the number of elements in Off , and d\ is 
the dimension of the A irreducible representation. Projections onto both the single-valued and double-valued irrcps 
of Oh require using the double group Off in Eq. (|36p . Given Mb elemental Bf operators, many of the projections in 
Eq. (|36p vanish or lead to linearly-dependent operators, so one must then choose suitable linear combinations of the 
projected operators to obtain a final set of independent baryon operators. Thus, in each symmetry channel, one ends 
up with a set of r operators given in terms of a linear superposition of the Mb elemental operators. The different 
spatial configurations (see Fig.[3]for the baryon configurations and Fig.Ufor the meson configurations) yield operators 
which effectively build up the necessary orbital and radial structures of the hadron excitations. The design of these 
operators is such that a large number of them can be evaluated very efficiently, and components in their construction 
can be used for both meson, baryon, and multi-hadron computations. 

Finding appropriate smearing parameters is a crucial initial part of any hadron spectrum calculation. Fig. [5] 
demonstrates that both quark-field and link-field smearing are needed in order for spatially-extended baryon operators 
to be useful^]. It is important to use the smeared links when smearing the quark field. Link smearing dramatically 
reduces the statistical errors in the correlators of the displaced operators, while quark-field smearing dramatically 
reduces the excited-state contamination. In this study, the quark field is Gaussian smeared with a = 3.0 using 32 
iterations, and the link field is stout-smeared with n p = 16,n p p = 2.5. 

Our approach to designing hadron and multi-hadron interpolating fields leads to a very large number of operators. It 
is not feasible to do spectrum computations using all of the operators so designed; for example, in the G\ g symmetry 
channel for nuclcons, the above procedure leads to 179 operators. It is necessary to 'prune down the number of 
operators. After much exploratory testing and trials, we found that a procedure that keeps a variety of operators 
while minimizing the effects of noise works best for facilitating the extraction of several excited states. Some operators 
are intrinsically noisy and must be removed. In addition, a set of operators, each with little intrinsic noise, can allow 
noise to creep in if they are not sufficiently independent of one another. 

The following procedure was used. (1) First, operators with excessive intrinsic noise were removed. This was done 
by examining the diagonal elements of the correlation matrix and discarding those operators whose self-correlators had 
relative errors above some threshold for a range of temporal separations. A low-statistics Monte Carlo computation 
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on a reasonably small lattice was used to accomplish this. (2) Second, pruning within operator types (single-site, 
singly-displaced, etc.) was done based on the condition number of the submatrices 

Ctf(i)= Vc«(»)' t = at ' 

The condition number was taken to be the ratio of the largest eigenvalue over the smallest eigenvalue. A value near 
unity is ideal. For each operator type, the set of about six operators which yielded the lowest condition number of the 
above submatrix was retained. (3) Lastly, pruning across all operator types was done based again on the condition 
number of the remaining submatrix as defined above. In this last step, the goal was to choose about 16 operators, 
keeping two or three of each type, such that a condition number reasonably close to unity was obtained. As long as 
a good variety of operators was retained, the resulting spectrum seemed to be fairly independent of the exact choice 
of operators at this stage. Eigenvectors from a variational study of the operators could also be used to fine tune the 
choice of operators. 



VII. 1=1 BARYON SPECTRA 

Each of our unbarred (barred) baryon operators annihilates (creates) a baryon and creates (annihilates) an an- 
tibaryon. This causes correlation functions to have a baryon state propagating forward in time and an antibaryon 
state propagating backward in time. Because fermions and antifermions have opposite intrinsic parities, the backward 
in time signal corresponds to states with parity opposite to that of the states propagating forward in time. Because 
of PCT symmetry coupled with the use of antipcriodic boundary conditions, correlation functions obey the rule 

c$® = c%;\T-ty, (a?) 

where A c is the opposite parity partner of irrep A. This allows us to increase our statistics by "folding" the correlation 
functions: 

- \ (C$(f) + C&°(T - t)*) . (38) 

Generally the separation of the two time slices involved is sufficient to provide independent samples of the gauge 
configurations. 

We choose phases of our baryon operators such that the matrices of correlation functions are real. We also average 
the matrix with its transpose in order to guarantee that the matrices are symmetric. This helps to clean up the 
signals by reducing the errors. 



A. The Variational Method 



We calculated 16x16 matrices of correlation function in each irrep of the octahedral group: A = 
{Gig, G2 g , H g , Giu, G2u, H u }. The variational method was used to help extract the excited spectrum from the matrices 
of correlation functions, which involved numerically solving the generalized eigenvalue problem 

C$(t)v$>(t,to) = ^n A) (t,t )Ciikto)v^(t,t ), (39) 

where n labels the eigenstates. Degeneracies and numerical uncertainties can cause variances of the eigenvectors 
at different times. We studied two methods for extracting the spectrum: 1) a fixed-eigenvector method and 2) a 
principal-correlator method where the diagonalizations were performed on each time step. 

The fixed-eigenvector method involved solving the eigenvalue problem on a single time slice t = t* using a fixed 

value of to- These eigenvectors are normalized with respect to C^J^to) such that 

v^(t*,t )Cit}(t )v^(t*,to) = S kk , (40) 

For each time slice and each configuration, the matrix of correlations functions was rotated to this basis of vectors 
using 

Cltkt) = V k l(t*,to)C^)(t)V n , k ,(t*,t ), (41) 
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where V n 'k'(t* ,to) is the matrix whose columns are the eigenvectors at time t*. Because of Eq. (|40|) . the rotated 
matrices of correlation functions are equal to the identity matrix at time t . 

The diagonal elements of the rotated correlation matrix are related to the energies by 

C$ (t) ~ e - B *<fi-*°) + a n e- E ^-^ + o(e^ E ^- E ^), (42) 

where the sum over terms involving a n vanishes at time t* , but can contribute away from t* . The term involving the 
first omitted energy, En+i, has been derived by Blossier ct al. (4p| . We extract the low lying energies by performing 
fully correlated ^-minimization fits, modelling the k th diagonal element of the rotated correlator matrix as 

= i 1 ~ A)e- Ek(t - to) + Ae- E ' {t - t0 \ (43) 

where Ek is the energy of the k th state. The second exponential captures the contribution of the higher energy states 
and allows us to fit the correlators to early time slices. The choice of coefneents in front of each exponential enforces 
Cl^{to) = 1, as guaranteed by Eq. (|4"0|) . We assume that the a n in Eq. [52] are negligible in the time range over which 
we perform the fit. 

We optimized our choice of to and t* using a method adapted from that in Ref. [171 ] . In that work, the optimal 
choice of to was determined for the extraction of the charmonium spectrum using a principal-correlator analysis. This 
optimal choice balanced the need for the contributions of higher energy states to have decayed away (suggesting larger 
values for to) and for the correlator to have a low level of noise (suggesting smaller values for to). The energies in the 
low lying spectrum were extracted by fitting the principal correlators for various values of to- For each value of to, 
the correlator was reconstructed from these fit energies and the eigenvectors using the spectral decomposition of the 
correlator matrix 

c tJ (t) = (0<(t)0 3 -(o)) = J2 ^^ ermat - ( 44 ) 

a 

The overlap factors Zf = (O|0j|a) arc related to the eigenvectors of the correlator by 

Zf = (V-^V^e" 1 ^ / 2 . (45) 

A % 2 -like quantity was defined to measure how well the reconstructed correlator described the original correlator 
matrix: 



x 2 



±N{N + l)(t max - t ) - ±N(N + 3) 



E E (CiA^-cirmc^MiCij^-cip-it')), m 



where (t,t') is the correlation matrix for the correlator CV,. Although the principal-correlator method actually 
yields time dependent overlap factors Z(t) (because the correlator matrix is diagonalized on all time slices), it was 
observed that the Z(t) were reasonably constant and the reconstruction was done using a single Z(tz) chosen at a 
time such that % 2 was minimized. For tz > to, the variation in % 2 as a function of tz was minimal. 

In this work, we adapt this technique for the fixed eigenvector method, finding optimal values for to and t*. We 
extract the 16 lowest energies in the spectrum by fitting the diagonal elements of the rotated correlator matrix, 
Eq. ([4Tj) . obtained using a range of values for to and t*. Reconstructing the correlator from these masses and the Z 
factors at t* , we choose the to and t* which minimize the x 2 ■ 

To correctly extract the energy spectrum, it is also crucial to select an appropriate range of time slices on which to 
fit the correlator. In particular, we would like to avoid time slices where the opposite-parity backward-propagating 
state contributes to the correlator. For mesons, where the forward and backward-propagating states have the same 
parity, the variational method simultaneously diagonalizes the forward and backward-propagating parts of a meson 
correlation function. This is not the case for baryons where the forward and backward-propagating states have 
opposite parities and different energies. The forward-in-time signals dominate at small values of time but they decay 
exponentially and the backward-propagating signals can become significant after some threshold value of time. We 
were able to extract the energies of the states by fitting the diagonal correlation functions using Eq. (|43[) without 
significant interference from the backward propagating signal for all channels except G\ u at m n = 416 MeV. In this 
channel, the backward propagating signal is dominated by the G\ g ground state, which is the lowest energy state in 
the spectrum. For our lattice at the lower pion mass, the backward-propagating G\ g signal decayed slowly enough 
and the temporal extent was small enough (due to the anisotropy) that the G\ u signals had significant backward 
contamination even at small time slices. To extract the G\ u energy levels using the fixed-eigenvector method, we 
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include the backward propagating state in the fit and constrain its energy by fitting simultaneously the G\ g ground 
state: 



= (i- A- B)e- E * lu + Ae- E '^-^ + Be E « 

C fit,Gi s = (1 _ D)e -E^(t-t ) + De -E"(t-t )_ 



'(t-t ) 



(47) 
(48) 



Due to the increased noise in the excited states, the minimizer was unable to find a minimum in the \ 2 for these 
simultaneous fits for k > 2. We were able to successfully fit these states by modeling the forward propagating state 
as single exponential and fitting only on later time slices (where the higher energy states had completely decayed). 

The fit ranges were optimized such that the x 2 was minimized. To visually confirm the sensibility of the fit 
parameters, we look at plots of 



^kk e > 



(49) 



versus time. If Eq. (|4"5|) correctly models the correlator, then the plot should plateau to (1 — A) and we confirm that 
the plateau is consistent with the value of A determined from the fit. For the G\ u channel at = 416 MeV, we 
first subtract off the backward exponential and compare the plateau with (1 — A — B) as in Eq. (|47|) . Finally, we 
confirm that the fit parameters are stable under small variations in the fit range. We estimate the uncertainty in the 
fit energy through a jackknife analysis. We fit each member of a jackknife ensemble to obtain an ensemble of energies 
and report the average energy and the jackknife error. 

The presence of the backward-propagating state in the G\ u channel caused numerical instabilities in the eigenvectors 
of the principal-correlator method. In order to remove the cause of the problem, we tested a method based on filtering 
out the backward signal prior to diagonalization. In a time interval where the backward signal is simply the ground 
state of the opposite parity channel with energy Eq" , the matrix of correlation functions can be modeled as a forward 
part plus the single backward state, 



-E'^(T-t 



(50) 



We define the filtered correlator as 



' filt,kk' 

and find that it can be modeled as 



(51) 



j=t+i 



C 



(A) 

filt.kk' 



A^ n ] = A^] 



l-e- E 'o c 
e E ~ - 1 



(52) 



where t\ is a time where the backward signal is, in fact, described by single exponential. The backward-in-time signal 
for energy E " is reduced to the level of errors and the filtered correlators consist of the renormalized forward signal 
minus a constant term. The diagonalization of the filtered correlators using the principal-correlator method produced 
stable eigenvectors and the energies of the states could be extracted by fitting the principal correlation functions to 
a single exponential decay with a constant term. However, this method did not produce any significant improvement 
over the results from the fixed-eigenvector method. We point out that the filtering is necessary in order to extract 
the Giu excited spectrum from our lattices using the principal-correlator method. 



B. Results 



We extracted spectra using the fixed-eigenvector method from the m„ = 416 MeV lattice using 430 gauge configu- 
rations and from the = 578 MeV lattice using 363 gauge configurations. Four states are reported for each channel 
for both pion masses. The results for = 416 MeV are given in Table |VT1 and the results for = 578 MeV are 
given in Table ["VTTl The results are based on 16x16 matrices of correlation functions using values of to, t* and the 
fitting windows U—tf as shown in the tables. Plots of the Nf = 2 spectrum for the two values are shown in Fig. 
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[H The pion mass is shown by the dashed line and thresholds for multiparticle states (to be discussed further on) are 
shown by empty boxes. Plots of Eq. (|49p versus time for each extracted state are shown in Figs. [7]fT2l 

In the positive parity channels, we identify the G\ g ground state as the nucleon. The spectrum for m OT = 416 McV 
is shifted toward higher energy values for = 578 MeV. The nucleon mass increases 172 MeV from 1136 MeV to 
1308 MeV when increases 162 MeV. If we extrapolate the nucleon to the physical pion mass using the formula 
M = a + bml, the result is 972(28) MeV. 

Results for the negative-parity excited states exhibit some interesting features. The pattern of G\ u energies shows 
two states at approximately 1.5 and 1.6 times the nucleon mass with the next state much higher. This pattern is 
similar to the pattern of masses of the physical spectrum, which has | resonances at 1535 MeV and 1650 MeV with 

the third \ resonance well above them at 2090 MeV. Because our baryon operators do not contain multi-hadron 
operators, they are expected to couple more strongly to three-quark states, suggesting that the lowest G\ u state is 
more likely to be a N* state. However, it is above the threshold for a ttN scattering state so further analysis clearly 
is needed to confirm this assignment. 

An isolated state in the H u irrep corresponds to a spin | state for which the lowest physical state is the iV(1520) 
resonance and the next to lowest is the 7V(1700). In the H u channel, the energies of the three lowest states are about 
1.57, 1.62 and 1.73 times the nucleon mass at the lower pion mass. The physical states for spin | are 1.62 and 
1.81 times the physical nucleon mass. In the G^u channel at to t = 416 MeV, we see that the lowest-energy state 
at 1957(51) MeV is degenerate (within errors) with the third H u state at 1964(48) MeV with no state at the same 
energy in the G\ u channel. A similar pattern is seen for = 578 MeV, except shifted upward by about 190 MeV. 
The lowest Gi u state at 2133(43) MeV is degenerate with the third H u state at 2182(38) McV. This pattern is the 
signature of a spin | state and a nearby spin | state. One H u state, most likely the second, is the spin | state 
and the other H u state is the partner state of the Gi u state required for a spin § state. The lowest possible spin 
in Gin is | and because the G^ u irrep has only two of the 2J+1=6 components needed for spin-|, the other four 
components necessarily are in a partner H u state. For a spin-| state, the Gi u and H u states must be degenerate in 
the continuum limit and for a clean interpretation there should not be a G\ u state that is degenerate with these two 
because that would be the signature of an isolated spin-| state or a possible accidental degeneracy of a spin i and 

| states. Our spectra show evidence for a spin-| state and a spin | state close to the same energy. As the pion 
mass is reduced to 140 MeV and the lattice spacing is extrapolated to zero, the partner H u and G2 U states in the 
lattice spectrum should approach the lowest | state in the physical spectrum, i.e., iV(1675) with a half- width of 75 

MeV. The first and second H u states should approach the 1520 MeV and 1700 MeV spin | states in the physical 
spectrum. 

The first excited positive-parity state in G\ g is at 2082(70) MeV for the lighter pion mass. That is 1.83 times the 
mass of the lowest G\ g state (nucleon) and about 334 MeV more massive than the lowest G\ u state. It also is well 
above the threshold energy for a p-wave Nir state(1785 MeV at the 416 MeV pion mass and this lattice length). In 

the physical spectrum the first excited, even-parity resonance is iV(1440)| with energy 1.53 times the nucleon mass 

and below that of the lowest odd-parity V*(1535)i state. Whether the energy of the first excited Gi g state will 
decrease toward the Roper state at lower values of the pion mass remains an open question. 

A signal for a | state could not be clearly identified in the quenched QCD analysis of Ref. @ at 480 MeV pion 
mass. That spectrum had larger errors and showed three degenerate states (within errors) in the Gi u , H u and G\ u 
irreps, a pattern with two possible interpretations. It could be a single spin-^ state or an accidental degeneracy of 

a spin-| state and a spin-i state. For Nf = 2 QCD and = 416 and 578 MeV, we see clear evidence for a | 
state. 

As the pion mass decreases, it becomes increasingly likely that some of the energy levels determined in our simula- 
tions will correspond to multi-hadron states. Disentangling these states from the hadron spectrum will be challenging 
and will require the use of specially-designed multi-hadron operators. In this paper, as a first step towards the iden- 
tification of scattering states, we estimate multi-hadron threshold energies in each of the irreducible representations 
of Oh- Some of the threshold energies correspond to states with two hadrons at rest. However, scattering states of 
hadrons with back to back momenta must also be considered. 

On the lattice, a hadron with momentum p transforms irreducibly under the space group, which is the semi-direct 
product of the group of three-dimensional lattice translations with Oh- In addition to the momentum vector p, irre- 
ducible representations of the lattice space group are characterized by a label denoting the irreducible representations 
of the group of lattice rotations which leaves p invariant (the little group of p). For particles at rest, the little group 
is Oh- More generally, the little group is a subgroup of Oh which depends on the orientation of p with respect to the 
lattice axes. The minimum non-zero momenta on a periodic lattice, of magnitude 27r/ (N s a s ), are directed along the 
lattice axes. The little group for such momenta is C\ v . 
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Gi 9 , io = 7, i = 10 Giu, to — 7, t — 9 

time £a t J5 (MeV) time Ea t E (MeV) 

3 - 21 0.2044(18) 1136(10) 3 - 14 0.3146(61) 1748(34) 

2 - 14 0.3747(126) 2082(70) 2 - 14 0.3343(67) 1857(37) 

2 - 12 0.4177(137) 2321(76) 7 - 14 0.5014(136) 2786(76) 

2 - 12 0.4201(277) 2334(154) 7 - 13 0.5238(158) 2910(88) 

H g , to = 8, t = 10 H u , to = 8, t = 9 

time Ea t E (MeV) time Ea t E (MeV) 

3 - 16 0.4004(74) 2225(41) 3 - 23 0.3208(87) 1782(48) 
3 - 17 0.4146(126) 2304(70) 3 - 21 0.3320(86) 1845(48) 
3 - 18 0.4193(120) 2330(67) 3 - 19 0.3535(87) 1964(48) 
3 - 16 0.4144(202) 2302(112) 2 - 11 0.5157(174) 2865(97) 
G2g, tp = 6, t* = 8 G2u, tp — 6, t* = 9 

time Ea t E (MeV) time Ea t E (MeV) 

2 - 12 0.4448(122) 2471(68) 2 - 17 0.3523(92) 1957(51) 

2 - 12 0.4593(104) 2552(58) 2 - 12 0.5035(119) 2797(66) 

2 - 11 0.4659(110) 2589(61) 2 - 12 0.5373(162) 2985(90) 

2 - 14 0.4796(127) 2665(71) 2 - 10 0.5446(131) 3026(73) 



TABLE VI: Isospin \ spectrum for m n =416 MeV. The energies in MeV units are based on the scale a t = 5556 MeV, and do 
not include the error in the the determination of the scale that acts as an overall multiplicative factor in the range 0.94 to 1.06. 





Gi 9 , to — 6, t 


= 10 




Giu, to — 6, t 


= 9 


time Ea t 


E (MeV) 


time Ea t 


E (MeV) 


2 - 


27 0.2463(17) 


1308(9) 


2 


- 11 0.3719(48) 


1975(25) 


2 - 


15 0.4291(110) 


2279(58) 


2 


- 11 0.3811(56) 


2024(30) 


2 - 


15 0.4643(116) 


2465(62) 


2 


- 11 0.5186(141) 


2754(75) 


2 - 


11 0.4631(123) 


2459(65) 


2 


- 11 0.5431(121) 


2884(64) 




H g , to = 6, t* 


= 9 




H u , to = 5, t 


= 7 


time Ea t 


E (MeV) 


time Ea t 


E (MeV) 


2 - 


- 14 0.4450(90) 


2363(48) 


2 


- 11 0.3802(86) 


2019(46) 


2 - 


- 11 0.4789(96) 


2543(51) 


2 


- 11 0.3975(89) 


2111(47) 


2 - 


- 11 0.4758(95) 


2526(50) 


2 


- 11 0.4110(72) 


2182(38) 


2 - 


- 11 0.4996(99) 


2653(53) 


2 


- 11 0.5670(215) 3011(114) 




G2g, to = 5, t 


' = 9 




Giu, to — 5, t 


= 9 



time Ea t E (MeV) time Ea t E (MeV) 

2 - 15 0.4422(144) 2348(76) 2 - 11 0.4017(81) 2133(43) 

2 - 15 0.4887(113) 2595(60) 2 - 11 0.5223(188) 2773(100) 

2 - 12 0.5030(94) 2671(50) 2 - 11 0.5399(139) 2867(74) 

2 - 14 0.5035(108) 2674(57) 2 - 11 0.5601(142) 2974(75) 



TABLE VII: Isospin \ spect rum for m v =578 MeV. The energies in MeV units are based on the scale a t 1 = 5310 MeV, and 
do not include the error in the scale determination that acts as an overall multiplicative factor in the range 0.95 to 1.05. 

Given the spectrum of hadrons at rest, one can deduce the allowed free-particle energies in any irreducible represen- 
tation of the space group. To see this, we first note that representations of a lattice little group can be subduced from 
the irreducible representations of Oh ■ The subduced representations are in general reducible and may be decomposed 
into a direct sum of irreducible little group representations. Irreducible representations of the full space group are 
induced from the irreducible representations of the lattice little groups. Thus, one can relate the irreducible repre- 
sentations of the space group to the representations of Oh- Neglecting cutoff effects, the energy of a non-interacting 
hadron with momentum p is given by E = yJ M? + \p\ 2 , where Mh is the rest mass of the hadron. Therefore, provided 
that the hadron rest masses are known, the free-particle energies in representations with non-zero p can be determined. 
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FIG. 6: The energies obtained for each symmetry channel of isospin | baryons are shown based on the 24 3 x 64 Nf = 2 lattice 
QCD data for m w = 416 MeV (left panel) and — 578 MeV (right panel). The scale shows energies in MeV and errors are 
indicated by the vertical size of the boxes. The overall error in the scale setting is not included. Empty boxes show thresholds 
for multi-hadron states. 
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FIG. 7: Plots of Eq. (|49|l versus time for Gi 9 states (left panel) and Gi« states (right panel) for m,r=416 MeV. For the two 
lowest energy states in the Gi u channel we first subtract off the backward exponential. 



Ref. [4l| gives the decomposition of direct products of irreducible representations of the space group, including 
representations with non-zero momentum, into the irreducible representations of Oh- We use this information to 
identify the allowed multi-hadron states in each representation of Oh- The energies of multi-hadron states are approx- 
imated by the sum of the energies of their constituents. The empty boxes in Fig. [6] show candidates for multi-hadron 
thresholds for both pion masses. Note that / = | baryons, which are not considered in this study, can also com- 
bine with isovector mesons to form / = | two-particle states. However, such states are expected to lie above the 
thresholds presented here. In both figures, the threshold energies in the G ltt and G 2 g representations correspond to 
meson-baryon states involving a pion at rest, while the other thresholds involve particles with non-zero momentum. 
The threshold energies in the Gi g , H g and H u representations arc degenerate. Our results illustrate the need for a 
proper analysis of multi-hadron contamination. Even at the heavier pion mass, many of the measured energy levels 
lie above the threshold for scattering states. Due to lattice artifacts, finite volume effects and the interaction between 
hadrons, the measured multi-hadron energies are expected to deviate from our estimates. This might explain some 
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FIG. 9: Plots of Eq. (|49[) versus time for G2 g states (left panel) and Giu states (right panel) for m^=416 MeV. 



of the discrepancies between the predicted multi-hadron energies and the measured spectrum. However, it is also 
likely that the interpolating operators used in our simulations, selected on the basis of a quenched study, couple only 
weakly to the lowest-lying multi-hadron states. Nevertheless, our analysis indicates that multi-hadron states cannot 
be discounted, even at the moderate pion masses used in this study. 
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VIII. CONCLUSION AND OUTLOOK 

In this work, anisotropic lattices with a t = |a s are developed for Nf = 2 QCD with two pion masses: = 416 
MeV and 578 MeV. The lattice setup and the algorithms used to generate gauge configurations arc described in detail. 
Conventional two-point correlation functions are used to calculate the spectrum of mesons in order to determine the 
pion masses and to tune the fermion anisotropy to £ = 3, which matches that of the gauge fields. The lattice scales 
a s ~ 0.113 fm and 0.108 fm are set using the Sommer parameter. 

This work builds upon several years of work to develop large numbers of baryon operators, to project them to the 
relevant irreducible representations of the octahedral group, to optimize the smearing of both the quark and gluon 
fields in the operators in order to be able to extract clean signals for effective masses and to prune the operators 
to manageable sets of 16 operators that yield good signals for baryons. Using the final operators, 16x16 matrices 
of correlation functions are calculated in each irrep and a variational analysis of the isospin | spectrum is carried 
out. The lowest four energy levels in each irrep are reported. The analysis of the negative-parity spectrum shows a 
cluster of states near 1.5 to 1.7 times the nucleon mass that includes a | state, two \ states at somewhat lower 
energies and two | states. This pattern is in accord with the pattern of physical states, although the latter is at 

a lower overall energy scale. The clear signal for a | state has not been realized previously. The analysis of the 
positive-parity spectrum for both pion masses shows that excited states typically have energies about 1.8 or more 
times the mass of the nucleon state. The question remains open whether as the pion mass is reduced the first excited 
G\g state will come down to about 1.53 times the nucleon mass, where it would agree with the Roper resonance. 

All the excited states in the lattice spectrum are near or above the threshold for irN scattering states. In order to 
deal properly with that aspect, multi-hadron operators and all-to-all propagators will be needed. This is an immediate 
challenge for progress on the 2+1-flavor dynamical lattices [lij [42| and it will be addressed in the near future. 
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